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We study the Hamiltonian dynamics of a free particle injected onto a chain containing a periodic 
array of harmonic oscillators in thermal equilibrium. The particle interacts locally with each oscil- 
lator, with an interaction that is linear in the oscillator coordinate and independent of the particle's 
position when it is within a finite interaction range. At long times the particle exhibits diffusive 
motion, with an ensemble averaged mean-squared displacement that is linear in time. The diffusion 
constant at high temperatures follows a power law D ~ y 5 / 2 f or a \\ parameter values studied. At 
low temperatures particle motion changes to a hopping process in which the particle is bound for 
considerable periods of time to a single oscillator before it is able to escape and explore the rest of 
the chain. A different power law, D ~ T 3//4 , emerges in this limit. A thermal distribution of parti- 
cles exhibits thermally activated diffusion at low temperatures as a result of classically self-trapped 
polaronic states. 

I. INTRODUCTION 

The theoretical problem of an essentially free particle interacting locally with vibrational modes of the medium 
through which it moves is relevant to a large class of physical systems, and arises in a variety of contexts. It is 
pertinent to fundamental studies aimed at understanding the emergence of dissipation in Hamiltonian systems 0, iM , 
and to the outstanding theoretical problem of deriving macroscopic transport laws from the underlying microscopic 
dynamics 0,0. It appears, perhaps most extensively, in the area of solid state physics [EIE0>I3j where a great deal 
of theoretical work has focused on the nature of electron-phonon interactions, and the rich variety of behavior that 
occurs: from the simple scattering of particles by phonons J5|, to the emergence of polarons in which the itinerant 
particle is clothed by a local distortion of the medium [y, [?], |8| , and even to self-trapped polaron states, in which the 
particle is bound to the potential well created by such a localized vibrational distortion. Such theoretical issues have 
also had a rich interplay with experiment, with experimental observations suggesting theoretically predicted band- 
to-hopping and adiabatic-nonadiabatic polaron transitions of injected charges carriers in organic molecular solids 
9j. 

Although much of the theoretical work on polaronic effects has appropriately focused on quantum mechanical 
calculations of polaron hopping rates, diffusion constants, and mobilities, there have been some workers who have 
questioned the degree to which polaronic phenomena depend intrinsically on their quantum mechanical underpinnings. 
There is an extensive literature, e.g., dealing with the validity and usefulness of various semi-classical approximations 
that have been used to study many aspects of polaron dynamics |lfj | . Considerably less work has focused on completely 
classical versions of the problem and the degree to which classical models mig ht capture essential transport features 
found in approximate or perturbative quantum mechanical calculations Moreover, since the exact quantum 

mechanical dynamics of the interacting particle/many-oscillator system remains computationally complex due to the 
enormous size of the many-body Hilbert space, it is difficult to computationally test perturbative or semi-classical 
calculations of such fundamental transport quantities as the polaron diffusion constant. It might be hoped, therefore, 
that numerically exact calculations of the full many-body dynamics of classical versions of the problem could provide 
information that would be complementary to analytical studies and semiclassical approximations of the full quantum 
evolution. 

Motivated by these considerations, we study here the closed Hamiltonian dynamics of a simple system that cor- 
responds to a classical version of the well-studied Holstein molecular crystal model 0. The model, in which a free 
mobile particle experiences a linear, local interaction with vibrational modes of fixed frequency lj arranged at regular 
intervals a along the line on which it moves, can also be viewed as a one-dimensional periodic inelastic Lorentz gas. 

Similar to what has been observed in quantum mechanical treatments of the polaron problem we find that there 
are in this model two different types of states of the combined system: self-trapped and itinerant. As their name 
suggests, self-trapped states are immobile. Itinerant states, by contrast, undergo motion that is macroscopically 
diffusive. We show for our classical model that the diffusion constant for itinerant states has a power law dependence 
on the temperature T, with two different regimes that are similar to the adiabatic and non-adiabatic limits discussed 
in the polaron literature. At high temperatures the diffusion constant of the model increases with temperature as 
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FIG. 1: Schematic illustration of the dynamical model. A particle moves along a line in which are periodically embedded 
oscillators of frequency ui which oscillate transverse to the particle motion. There is no interaction when the particle is in one 
of the free regions of length L between adjacent oscillators. The particle exerts a constant force on the oscillator at q m = ma 
when it is within a distance a of that oscillator. As a result, the oscillator in that cell moves about a shifted equilibrium 
position <j> t = —ol/uj 2 . The thick line shows a snapshot of the instantaneous potential U(q) experienced by the particle at one 
instant during the evolution of the system. The shaded regions connect the potential in each interacting subcell to the value 
about which it oscillates. 



T 5 / 2 and the transport mechanism is similar to what one thinks of as band transport in a solid, with long excursions 
before a velocity reversing scattering from one of the oscillators in the system. At low temperatures transport occurs 
via hopping of the carrier between different unit cells, with the particle spending a long time in the neighborhood of 
a given oscillator before making a transition to a neighboring one. In this hopping limit the diffusion constant for 
untrapped particles varies as T 3 / 4 , ultimately vanishing at zero temperature. An ensemble of particles in thermal 
equilibrium with the chain contains both itnerant and self-trapped populations of particles, and exhibits a diffusion 
constant D ~ T 1 ^ 2 e~ T °^ T that is exponentially activated at low temperatures, with a functional form similar to that 
derived for polaron hopping in quantum mechanical treatments of the problem. 

In the current classical model, we take the particle-oscillator interaction to be linear in the oscillator coordinate 
when the particle finds itself within a distance a < a/2 of the oscillator. Thus, the line on which the particle moves 
naturally decomposes into unit cells of width a, each centered on an interacting subcell of width 2a. This interacting 
subcell is surrounded on each side by noninteracting regions of width a — 2a = L. (See Fig.l.) The total Hamiltonian 
describing this system can be written 

H = ^P 2 + \ (^m + u *<t>m) + a Yl (<?) (1) 

tci m 

where p and q are the momentum and position of the particle, and n rn and <p m the momentum and displacement of 
the m th oscillator, which is located on the chain at q — q m — ma. In the Hamiltonian Q), the function n m (q) = 
8 (a — \q — q m \) governs the range of interaction between the particle and the mth oscillator: it vanishes outside the 
interaction region associated with that oscillator and is equal to unity inside it. We denote by 9 (x) the unit increasing 
step function. The parameter a governs the strength of the interaction. 

Although not necessary, it is convenient to think of the oscillator motion in this model as occurring in a direction 
perpendicular to that of the particle, as depicted in Fig.^ At any given instant the interaction term defines a 
disordered step potential U(q) = a^ TO (p m n m (q) seen by the particle as it moves along the chain. Under the classical 
evolution of the model, the vibrational modes in subcells in which the particle is absent (i.e., modes for which 
n m (?) = 0) oscillate independently around their non-interacting equilibrium position cp m = 0. During times that the 
particle is in one of the interacting regions, however, the associated vibrational mode in that region oscillates about 
a displaced equilibrium position 0* = ~a/uj 2 . This local displacement of the equilibrium position of the oscillator 
in the region around the particle can be viewed as a vibrational distortion of the medium that "accompanies" the 
particle as it moves, forming part of this classical polaron. Indeed, the ground state of this model corresponds to a 
polaronic "self-trapped" state in which the particle is at rest in one of the interacting regions, with all oscillators at 
rest in their normal or displaced equilibrium positions. For this state, the total energy of the system 
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is negative and has a magnitude Eb = a 2 /2u! 2 , analogous to what is referred to in the polaron literature as the 
"polaron binding energy". Thus, the binding energy Eb provides a natural measure of the coupling between the 
particle and the local oscillator modes. 

In our classical model, except for instances when the particle is at the boundary between an interacting and non- 
interacting region, the particle and oscillator motion are decoupled and evolve independently, with the particle moving 
at constant speed. At "impacts" , i.e., at moments when the particle arrives at the boundary between a non- interacting 
region and an interacting region, the particle encounters a discontinuity AU = A = ±a(j) m in the potential energy, 
and thus an impulsive force, that is determined by the oscillator displacement in the interacting region at the time 
of impact. When the particle energy e at the moment of impact is less than A, the particle is reflected back into the 
subcell in which it was initially moving; when e > A the particle moves into the next subcell, with a new speed and 
energy Ef = e — A governed by energy conservation. Thus, the general motion of the particle along the oscillator 
chain involves sojourns within each interacting or noninteracting subcell it visits, during which it can reflect multiple 
times, interspersed with impact events in which the particle's energy and the oscillator displacements allow it to pass 
into a neighboring cell. 

This uncoupled evolution of the two subsystems between impacts makes it easy to perform a numerical "event 
driven" evolution of the system by focusing on the changes that occur during impacts. We have performed a large 
number of such numerical evolutions of this model for a set of initial conditions in which the particle is injected 
into the chain at its midpoint, with each oscillator initialized to a state drawn independently and randomly from 
a thermal distribution at temperature T. Our numerical studies, which are supported by theoretical analysis given 
below, suggest that for such initial states there are three distinct types of dynamical behavior that can occur. One 
of these is essentially trivial, and leads to localized self-trapped particle motion only. The other two types both give 
rise to diffusive motion of the particle over large time and length scales. In all cases the oscillator distribution as a 
whole remains in thermal equilibrium. 

The self-trapped states of the model correspond to configurations in which a particle is located initially inside an 
interaction region (at q m , say), with the energy 

Emt = \p 2 + \ (^m + ^m) + a §™ ( 3 ) 

of the particle and the oscillator in that region partitioned in such a way that the particle is unable to escape. Such 
states include all those bound states in which eint is negative. As shown in Ref. however, it also includes a set 
of positive energy self-trapped states, i.e., all dynamical states for which Eb > £int > 0, and for which the associated 
oscillator energy 

£osc = \ (f m + o? (0 m - <E c = E B -e(\- -^A (4) 

is sufficiently low that the potential barrier seen by the particle at each impact is always greater than its kinetic energy. 
Clearly for all such initial conditions the mean square displacement (q 2 (t)) = (rn 2 (t))a 2 of the particle, measured (or 
course-grained) in terms of the mean square number (m 2 (t)) of unit cells visited, remains equal to zero. 

Our numerical studies indicate that for initial states in which the particle is not in a self-trapped state, the particle at 
sufficiently long times undergoes diffusive motion, with a mean square displacement (q 2 (t)) ~ Dt that grows linearly 
in time. Figure [21 shows a representative set of numerical data for the (dimensionlcss) mean square displacement 
(m 2 (t)) = (q 2 (t))/a 2 for several different parameters and temperatures T, the latter expressed in terms of the 
reciprocal temperature (3 — 1/kT, where k is Boltzmann's constant, and with each curve corresponding to particle 
evolutions (or trajectories) averaged over N t = I0 5 independent initial conditions |13|. In the numerical calculations 
presented here, time is expressed in terms of the oscillator frequency through the dimensionlcss combination r = cut, 
lengths are measured in terms of the width 2a of the interaction region, and energies are measured in terms of the 
quantity 

E Q = a 2 u 2 , (5) 

which is one-half the kinetic energy of a particle that traverses an interaction region in time t = uj^ 1 . Thus, aside from 
the temperature, there are only two independent mechanical parameters in the problem, which we can take to be the 
ratio r] — 2a I L of the width of the interacting and noninteracting regions, and the ratio e = Eb/Eq = a 2 /2a 2 io A of 
the polaron binding energy to Eq- 

The diffusion constant D describing the linear increase in the mean square displacement of initially untrapped 
particles, obtained from a fit to the numerical data for (q 2 (t)) in the region where it becomes linear, is found to be a 
strong function of the mechanical parameters and of the temperature. Results for a range of parameters and a wide 
range of temperature are presented in Fig. [3] Note that data for different values of the coupling strength parameter 



4 




FIG. 2: Mean square displacement (measured in units of the squared cell spacing a 2 ) for an ensemble of particles injected 
with positive energy into an oscillator chain in thermal equilibrium at temperature T = l/k/3, with temperatures, coupling 
strengths, and mechanical parameters as indicated. For the parameters in the linear plot on the left (m 2 (t)) appears linear 
over the entire range of times. On the log-log plot on the right one can see that linear growth of (m 2 (t)) begins after an initial 
period of super- or sub- diffusive behavior. The temperature range in the latter figure covers both the high temperature and 
low temperature regimes analyzed in the text. 

Eg I Eq are shifted by powers of ten to separate the data. If this scaling were not done the data would all lie on a single 
pair of universal curves. As the figure shows, for any given fixed set of parameters, the temperature dependence of the 
diffusion constant takes two simple limiting forms. At high temperatures, the dependence of the diffusion constant 
on temperature follows a power law 

D/D% ~ (fJE B )- 5/2 (6) 

where D° H is a constant that depends on the mechanical parameters. In Sec. |^we derive an approximate expression 
for D° H [See Eq. 1)21(1] which we have used to scale the data on the left hand side of the figure. The straight lines in 

the left hand side of the figure are the function {(3Eb) 5 ^ 2 ■ 

Also as seen in Fig. [21 as the temperature is lowered the diffusion constant crosses over to a weaker power law of 
the form 

D/D° L ~ (f3E B y 3/i (7) 

where D° L is a different constant, which we derive in Sec. IIIII below and which we have used to scale the data on the 

right hand side of that figure. The straight lines in the right hand side of that figure are the function (0Eb)~ 3 ■ 

These two different power laws result from two quite different forms of microscopic dynamics that, in a sense, 
correspond to the adiabatic and nonadiabatic polaron dynamics studied in quantum mechanical treatments of the 
problem. Indeed, at sufficiently high temperatures the characteristic particle speed p ~ s/kT becomes sufficiently 
large that the time for a particle to traverse the interaction region is short compared to the oscillator period. In this 
limit, also, the typical potential energy barrier A ~ \/2EskT encountered by the particle at impacts becomes small 
compared to the particle energy. The picture that emerges is that, at high temperatures, a particle with a (high) 
initial speed, close to the thermal average, will pass through many interaction regions in succession, typically losing 
a small amount of energy to each oscillator as it does so. Eventually, this continual loss of energy, which acts like 
a source of friction, slows the particle down to a speed where it is more likely to encounter a barrier larger than 
its energy, at which point it undergoes a velocity reversing (or randomizing) kick back up to thermal velocities. By 
heuristically taking into account a typical excursion I (p) of a particle with initial momentum p we derive in Sec. [n] 
the result © seen in our numerical data, and in the process derive an approximate expression for the scaling prefactor 
Djj. This regime is similar to the so-called adiabatic limit discussed in the small polaron literature, in which the time 
for a bare particle to move from one cell to the next is short compared to the evolution time of the oscillator. Indeed, 
in this limit the random potential seen by the particle typically changes adiabatically with respect to the particle's 
net motion. 

As the temperature is lowered, however, the mean particle speed decreases. At low enough temperatures a typical 
particle traverses the interaction region in a time that is large compared to the oscillator period, so that at each 
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FIG. 3: Diffusion constant obtained from linear fits to long-time mean square displacement, for a wide range of system 
parameters and coupling strengths as indicated, as a function of the inverse temperature (5Eb- Curves corresponding to 
different coupling strengths are multiplied by different powers of 10 as indicated for clarity. Data on the left (right) hand side 
has been scaled by the high (low) temperature value D° H {D° L ) derived in the text. Straight lines are the power laws (/3Eb) 5 ^ 2 
for the high temperature data and (/3i?s) 3 ^ 4 for the low temperature data. 

impact the oscillator phase is randomized. In this limit, also, the typical barrier energy \J2EskT becomes large 
compared to the particle energy. As a result, a particle will typically undergo many reflections inside the interaction 
region before escaping into one of the adjacent regions. The resulting motion can be described as a random walk 
between neighboring unit cells, with the diffusion constant depending on the escape rate from a given cell. In this 
low-temperature hopping limit, the motion is similar to the nonadiabatic limit of small polaron motion, in which the 
time to move freely from one cell to another is very long compared to the typical oscillator period. In Sec. 11111 we 
calculate this escape rate in order to produce a theoretical estimate of the diffusion constant in the low temperature 
nonadiabatic limit, represented by the solid lines on the right hand side of Fig. [3J 

In Sec. IIVI we combine the results of the two analyses in the preceding sections to show that the behavior of an 
ensemble of particles in thermodynamic equilibrium on such a chain would exhibit diffusive behavior that retains 
all essential features of the high temperature limit, but becomes exponentially suppressed in the low temperature 
hopping limit as a result of the increasing fraction of self-trapped particles at low temperatures. This results in a 
thermally activated diffusion constant very similar to that derived in the nonadiabatic limit using quantum mechanical 
arguments. The last section contains a summary and discussion. 

II. DIFFUSION OF INJECTED CARRIERS AT HIGH TEMPERATURES 

In this section we consider the dynamical behavior of this model at high temperatures, i.e., the left hand side of 
Fig. © . We note first that when a particle which is initially in one of the non-interacting regions attempts to enter 
an interacting region, the distribution of potential energy steps A = a<p m it encounters is given by the relation 

"< a > = S" a "' 4& - <8 » 

which follows from the assumed thermal distribution p((f> m ) = \J (3ui 2 /2ire~P ul 'I'm./ 2 for the displacement cf> m of the 
oscillator at the moment of impact and the definition J2J. Consequently, the magnitude of the typical (RMS) potential 
step 

A = ^2E B kT (9) 

seen by the particle, independent of its energy, increases only as the square root of the temperature. At asymptotically 
high temperatures the typical step is, therefore, relatively small compared to the energy e ~ kT/2 of the particle. 
Thermalized particles will therefore easily pass through many cells in the same direction with high probability. In 
this limit the traversal time becomes short compared to the oscillator period, and the particle will tend to lose energy 
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to the oscillator bath at a well defined average rate. To see this, we consider a particle with an initially high initial 
momentum p which has just entered the interaction region associated with one of the oscillators in the chain. Before 
impact, the oscillator in that region is (by assumption) in thermal equilibrium about the undisplaced equilibrium 
position (4> — 0)j an d so the average energy step (Ao) = ot(4>a) encountered by the particle at the interface vanishes. 
After entry, the state of the oscillator in the interaction region is characterized by the distribution 



P(^o) -sj^^l^l-e-^e (p 2 - 2c^ ) , (10) 

where the upper cutoff on <fi associated with the step function arises from the requirement that the initial energy of 
the particle and the oscillator state be such as to have allowed the particle to enter the interaction region in the first 
place. For high temperatures (i.e., 0E B <C 1) and particles with thermal energies p 2 ~ kT, this cutoff 4> c ~ kT/2a is 
well outside the range <pth ~ VkT/u of the thermal distribution for the oscillator displacement and can be neglected, 
i.e., <p c /4>th ~ {(3E B ) 1 3> 1. Once the particle is in the interaction region, the oscillator state, and hence the 
distribution p (<j), n, t) will evolve in time as the oscillator starts to oscillate about its displaced equilibrium position 
0*. After a time t = 2a/p, the particle will have traversed the interaction region, and the oscillator coordinate <j> will 
have changed by an amount 

Acp = (0o - </>*) (cos(2wcr/p) - 1) + — sin(2wer/p) . (11) 

D 

Averaging this over the initial distribution (110ft gives the average change in the oscillator displacement 

{A<f>) = 0, (1 - cos (2ua/p)) (12) 

at the moment the particle attempts to leave the interaction region after one pass. When it does so, the average 
difference between the energy steps that the particle crosses between entering and leaving the interaction region must 
equal the average change (Ae) in the particle's energy as it passes through that region, i.e., 

(Ae) = -a{A(f>) = -2E B [1 - cos (2ua/p)] . (13) 



Of course, at temperatures high enough that VkT ^> 2au, the oscillator period will be large compared to the 
characteristic time for particles of average thermal energy to traverse the interaction region. For such particles, we 
can expand the cosine in l|13fl to find that the average energy change 

(Ae) = -AE B E /p 2 (14) 

is negative, simply because the particle exerts a constant force on the oscillator in the interaction region through 
which it is passes that tends to decrease the potential energy in that region while it does so, and as a result makes 
it more likely for the particle to lose net kinetic energy when it passes back out into the noninteracting region. Thus 
the particle tends to slow down in time as it passes through many such cells in succession. The average momentum 
loss 

dp _ dp (Ae) _ 4E B E ^ 
dn de dn p 3 ' 

per cell traversed can then be integrated to obtain the typical distance 

t(p)=an(p) = -^— (16) 
16E B E 

the particle must travel, losing energy at this rate, to dissipate most of its initial momentum p. This is, in fact, an 
overestimate since the particle need only dissipate as much momentum is required for it to encounter, with significant 
probability, a barrier larger than its (reduced) energy. Of course, in the high temperature limit, when barriers of 
sufficient energy are few, the particle will need to dissipate a significant fraction of its initial energy. We also note 
from l|15|) that for fast enough particles the fractional change in momentum per traversal is small, so the average time 
to traverse each complete cell can be written dt/dn = a /p. This allows us to integrate the average "frictional force" 

dp _ dp^dn _ 4E B E ^ 
dt dn dt p 2 a 
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on particles with large initial momentum to obtain the characteristic time 

^=mk (18) 

for this slowing down to occur. 

The picture that develops from this analysis, therefore, is that, fluctuations aside, at high temperature a thermal 
distribution of particles will contain many particles of high energy that are slowing down on average in response to a 
frictional force arising from their interaction with the oscillators through which they pass. This slowing down can not 
continue, however, or else the distribution would not maintain itself in thermal equilibrium. Indeed, once the particle 
speed is slow enough that it becomes likely to encounter a barrier larger than its kinetic energy it receives a velocity 
reversing impulse. After such an event, if the time for the particle to traverse the interaction region is still small 
compared to the oscillator period the particle will continue to lose energy, on average, until its traversal time becomes 
comparable to the oscillator period, in which limit it can gain as well as lose energy due to the random phases of the 
oscillator at the moment it attempts to leave the interaction region. In the overall process a particle with a given 
initial momentum p will travel a distance £ (p) in a time r (p) before ultimately being kicked back up to higher energy. 
The associated diffusion constant can, we argue, be calculated as in a random walk with a distribution of hopping 
times r and step lengths £, i.e., 

f 00 £ 2 (v) 

D = 2 dpp(p) 1 -^- (19) 



r(p) 



where p (p) — /3/27te~ l3p I 2 is the thermal distribution of particle speeds. Assuming that the main contribution to 
the growth of the mean square displacement comes from the fast particles with long excursion lengths, as described 
above, we combine Eqs. I|18|) . and l|19|) to write the diffusion constant in the high temperature limit as 

D = 3a «/A r p 5 e~^ 2 / 2 dp = D° H ■ (f3E B y 5/2 (20) 
32E B E \ 2nJ y H \H B ) V ) 

where 



o 9E B a 2 E B , s 

This expression for the diffusion constant is represented as the solid lines in the left hand side of Fig. |3| Agreement 
of the numerical data with the scaling behavior predicted above is excellent. This collapse of the data for a wide range 
of values of E B /Eo and 2a/ L suggests that the heuristic arguments given above capture the essential physics taking 
place in the high temperature limit, even though they lead to an estimate of D Q H that is too low by a numerical factor 
lying between 1 and 2.5 over the range of parameters studied. 



III. DIFFUSION OF INJECTED CARRIERS AT LOW TEMPERATURES 



In this section we provide an analysis aimed at understanding the diffusion constant in the system at very low 
temperatures, i.e., the right hand side of Fig. J3J). As we have already argued, in this limit the barrier seen by 
a particle in a given subcell can be large compared to typical particle energies. Consequently, an injected particle 
that makes its way into an interaction region can rattle around in that region, undergoing many reflections before 
encountering a barrier sufficiently low that it is able to pass back out. Also, as the temperature is lowered, the typical 
time for a thermal particle to traverse the interaction region becomes very large compared to the oscillator period. 
Hence, after each traversal, the phase of the oscillator is randomized, except for a set of velocities of zero measure 
for which the traversal time is commensurate with the oscillator period. Thus, the motion of the particle at very low 
temperatures becomes a random walk between neighboring cells, and calculating the diffusion constant 

D ~ Ya 2 (22) 

reduces to a calculation of the average total escape rate T for leaving an individual unit cell. 

It is relatively easy to see, also, that at low enough temperatures the time for the particle to leave a given unit cell 
is dominated by the time for leaving the interaction region at its center. Intuitively this is clear, since as a thermal 
particle approaches an interface from the non-interacting side, the oscillator displacement in the interaction region 
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is equally likely to be positive or negative, but with a magnitude A ~ \JEskT that is large relative to the particle 
energy e ~ kT/2. Hence, at low temperatures the probability for the particle to pass into the interaction region 
at any given impact approaches one-half. Once inside the interaction region, however, the oscillator therein starts 
oscillating about its displaced equilibrium position. The initial amplitude of the oscillator at the moment the particle 
enters will typically be negative (so the particle could get in) and with a magnitude \4>o\ <~ ^/kT/to that is small 
compared to the magnitude of the maximum displacement max ~ —20* + |0o| ~ —2<ft« of the oscillations that it will 
now perform about its suddenly shifted equilibrium position. If, after the particle traverses the interaction region, 
the phase of the oscillator has become randomized, it is then very unlikely that the amplitude will allow the particle 
to exit the interaction region, since for most of the oscillation the amplitude is a significant fraction of ma xj and 
thus will present a barrier A which is a significant fraction of the maximum barrier A max ~ 2a<j)* ~ 4_Eg, which in 
the low temperature limit /3Eb 3> 1 is large relative to the particle energy. Consequently, in this limit the particle 
will typically reflect many times inside the interaction region until it arrives at the interface at a moment when the 
oscillator phase is very close to its value at the moment it entered. 

Calculation of the diffusion constant in the low temperature limit, therefore, reduces to a calculation of the average 
escape rate for a particle which has made its way into the interaction region from the outside. To this end, we 
focus on that subensemble of particles which have just entered an interaction region, and which end up with energy 
e = p 2 /2 once inside the region. As we will see, the conditional distribution P((fi,ir\ e) of oscillator states at the 
moment a particle enters the interaction region with this energy is straightforward to calculate, and allows us to 
obtain the corresponding conditional distribution P(e osc |e) of oscillator energies. Assuming that the phase of the 
oscillator is completely randomized as the particle traverses the interaction region, we can then use this to compute 
the distribution p (A| e) of potential energy barriers A = acp that a particle with this energy will see as it attempts to 
leave the region. It will do so if the barrier A is less than the energy e of the particle, so the a priori probability Q (e) 
that a particle which has entered the interaction region with energy e will exit the region after any given impact is 

Q (e) = / Pose (A | e) dA. (23) 

J — OO 

The average rate at which a particle of this energy leaves the interaction region can then be written as the product 

l(e) = Q(e)£ (24) 

of the attempt frequency v = p (e) /2a and the success probability Q (e). The total average rate T at which particles 
leave the interaction region can finally be evaluated 

T= T P i : t {e) 1 (e)de (25) 
Jo 

using the the steady state distribution p 1 ^ (e) of untrapped particles in the interaction region. 

We focus initially, therefore, on particles just about to ent er the int eraction region, with external initial energy e out , 
which is assumed to be thermally distributed, p{e out ) — \J /3/7re out exp (— /3e ou t ) , as is the state of the oscillator in 
the interaction region it is attempting to enter. Under these circumstances, the particle can only enter the interaction 
region if e out is larger than the potential step it encounters, i.e., if e out — a<t> = e is positive, where e is the energy of the 
particle once it has entered the interaction region. Up to a normalization constant, therefore, the joint distribution 
of the oscillator state and the particle energy at the moment it enters the interaction region is 



f(n,J>,e) = J g &e-WW)Pe-K***Me + <*l>) (26) 
y 7r (e + a<j>) 2ir 

where the step function ensures that the original particle energy and oscillator configuration were such as to allow it 
to enter in the first place. The corresponding conditional distribution of oscillator states in an interaction region into 
which a particle of energy e has just passed is then 



<> = ^ J Vra) e " ( " v+ ** ,/2e "" + "' S(E+a * ) (27) 



where 



Z(e) = fjcf, |~ dn f (tt, «M = ^1 - ^e-^- E ^e- 2 / 2 K 1/4 (« 2 /2) 



(28) 
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In this last expression, K v (z) is the modified Bessel function, and we have introduced the quantity 

e 



k (e) = V0E^ 1 



2E B 



(29) 



With the particle in the interaction region, the total mechanical energy of the oscillator as it oscillates about its 
displaced equilibrium position can be written 



l r 



l r 

2 



2 i 2 ±1 
7T + UJ <t> 



= - 7T 2 + ui 2 (4> ~ (j)*)' 

The conditional distribution for the oscillator energy can be computed from (|27|l using the relation 

1 



P(e osc \ e) 
We find, after some work, that 

P{e OS c\ e) = 



oo poo 



oo J — oo 



dn P (tt, <j)\e)6[ s 



A(e OBC ,e)K y7t( 

^OSC ) & ) 



2 I 2 72 
TT + UJ 



e < E- 

£_<£<£+ 



(30) 



(31) 



(32) 



yjfi (e osc ,e)A(e osc ,e) K ^ n" 1 {e osc ,e) 



e + < e 



where K (m) is the complete elliptic integral of the first kind, e± = 2 (Eb ± \/e osc Eb) , and we have introduced the 
functions 



A (s osc , e) 



^3/2 



tt^ 2 {e osc E B ) l/i Z (e) 



-/3(e+e oac -E B ) 



and 



A 4 (£o SC ,e) 



Ve osc /E B + e/2gg - 1 
2 \/ £osc / 'Eb 



(33) 



(34) 



As we will see below, at very low temperatures the distribution p 1 ^ (e) of untrapped particles in Eq. H25fl . becomes 
exponentially suppressed at higher energies, and it suffices to consider the form of P (e OS c|e) for values of e less than or 
of the order of a few kT <C Eb- For such values of e, the function P (e osc \ e) , as a function of e OS c, becomes strongly 
peaked in the neighborhood of e osc ~ Eb- This can be seen from Ea. (|32|l . where it is clear that for values of e in this 
range the condition for validity of the last line is never satisfied. The first line then shows that P(e OS c| e) vanishes 
except for oscillator energies e osc satisfying 



> E B - e 1 - 



4E f 



E B [1- 



1 



(3E E 



(35) 



which is very close to, and just below, Eb at low temperatures. Moreover the function A (e osc , e) drives the distribution 
P (£osc| e) exponentially towards zero for increasing oscillator energies greater than a few kT of where it first becomes 
non-zero. Thus, at low temperatures, for values of e of interest, the distribution P(e sc| s) becomes increasingly 
peaked in the variable e osc in the neighborhood of Eb- Since the function /i (e sc,£) ~ s/4Eb is small in this region, 
we can approximate the elliptic integral in (1-i2l) by its limiting value K (0) = n/2. After some analysis, we find that 
at low temperatures the distribution P(e OS c| £) asymptotically approaches the limiting distribution 



P (e OS c| e) 



-pe osc 



£ — 2 (Eb — \/ £ osc Eb 



(36) 



where 



„_r(3/4, K 2 ) 



2\ ex P 



-0E B (l 



2E B 



/33/ 4 / , N 2\V4' 

^ ^(1-31-, 



(3E B > 1, 



(37) 
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in which T (a, z) is the incomplete gamma function, the asymptotic form of which we have used in the second expression. 

Now for a given oscillator energy £ OS c> an d assuming uniformly distributed phases, the conditional distribution for 
the displacement <f> of the oscillator about its displaced equilibrium position 0* , given the oscillator energy e osc is easily 
computed. Using this, the distribution of barriers A = — acf) encountered by a particle after traversing an interaction 
region in which there is an oscillator of energy e osc is readily found to be 



p B (A| e Q 



+ 0*)' 



2tt 



E b (e 



j 2 l A 



(38) 



Multiplying i|36|) by i|38|l an d integrating over all oscillator energies allows us to obtain the conditional barrier 
distribution 



p(A\e) 



cfeosc P (e OSC \ e ) Pb (A| 

c osc ) 



where, due to the step functions in l)36|) and (|38[) . the lower limit of integration 

E\ — XEb 



A = max 



1 



2E B 



A 
2Eb~ 



1 



(39) 

(40) 
(41) 



in i|39|) is a function of e and A. As we have already observed, at low temperature the integrand in (|39|) dies away 
exponentially as a function of e osc from the maximum value it takes on at the lower limit. By expanding to lowest 
order the non-exponential part of the integrand about the lower limit of integration and performing the resulting 
integral we find that for j3Es ^> 1 and for values of e < 2Eb, the function p(A\ e) is well approximated by the 
expression 



P(A| e) = 



2nE E 



1 2E B J \2E B 1 



-1/2 







V2E B -e „-PEt 



-I3E E 



4ttE b \A-2E b \ 1 ' 2 

From the barrier distribution 1421) we can then compute the probability 



for 
for 



1_^ 2 > P^-l 

1 2E B J ^ \2E B 1 



(42) 



A 

2E B 



> 1 



2E E 



•(e) 



dA P (A\e) 



(43) 



that a particle that has entered the interaction region with energy e will pass out of it as the particle approaches 
the interface. Since the integrand in (|43|) only includes values of A less than or equal to e, only that part of the 
distribution 142(1 on the second line of the expression is needed. Using this in (|43(l the resulting integral gives for 
e,P~ 1 <$:2E B 



'(e) 



T(l^)(f3E B ) 



1/4 



4tt 



1 - 



-exp 



1 



2E B 

(3E B > 1, E B ^> e 



&Eb 1 - 



2E B 



(44) 



where in the last line we have again used the asymptotic expansion of the incomplete gamma function for large 
arguments. 

The last distribution needed to calculate the diffusion constant using Eqs. (1221 and (|25|l is the steady-state distri- 
bution p™ (e) of untrapped particles inside the interaction region. Starting with the canonical distribution 



P mt (e OS c,e) 



e ^ s e-^°- 



(45) 



of the combined particle-oscillator system inside an interaction region, we exclude that portion associated with trapped 
states by multiplying by an appropriate step function that incorporates the self-trapping condition of Eq. i.e., 



//J3 

(^osci^) — \/ *2 
V 7T£ 



■Ps e -l3e os 



-E B +e 



1 - 



4E B 



(46) 
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We now integrate over the oscillator energies to find the associated particle distribution 



± e -PE Be -^/4E B 2E B >e>0 
C (e) = PT (e OS c, e) ds osc = {V™ (47) 

-£- e -" B £ > 2£ B 



which is normalized to the fraction /™* of untrapped particles in the interaction region in thermal equilibrium. At 
very low temperatures (3E B 3> 1, the fraction of untrapped particles with energy e > 2E B becomes negligible and we 
can approximate the normalized distribution p u (e) of untrapped particles in the interaction region for all positive e 
by the form it takes for e < 2E Bl i.e., 

i„t _ /r (f) r(3/4) / g \ 1/4 » £ * /4Ee 



PT 00 = f oc7 int V w — -r^r- e"^ . (48) 

Jo fuH^de 7T \e 2 E B J 

Using Eqs. (|24|) . I|25|) . I|44|) . and l|48|) the average escape rate can be written 

p = — f°° V2lo int (e) O (e) de - F [°° e -^l^ B de 



Using the resulting diffusion constant then takes the form 

D = D° L (pE B )- 3/i (49) 

where 



^ = i; r < 3 / 4 >V&- 



As seen in Fig. [21 this functional form does an excellent job of describing the temperature dependence of the diffusion 
constant for injected charges at low temperatures for a wide range of model parameters. 



IV. DIFFUSION OF CARRIERS IN THERMAL EQUILIBRIUM 



The previous sections focus on calculating the diffusion constant for particles injected into an untrapped state in a 
thermalized chain of oscillators. An ensemble of noninteracting carriers in thermal equilibrium with the oscillators in 
the chain will contain both itinerant particles, which will exchange energy with the oscillator system as they undergo 
diffusion along the chain, as well as self-trapped particles, with energies satisfying Q and which remain bound to the 
same oscillator forever. Since the diffusion constant for untrapped particles is identically zero, the ensemble averaged 
mean-square displacement of a collection of particles in thermal equilibrium will be characterized by a diffusion 
constant that can be written as the product 

D = D u f u (51) 

of the diffusion constant for untrapped particles, i.e., the diffusion constant evaluated in the last two sections for high 
and low temperatures, and the fraction 

/>oo 

fu = / fu (e) de (52) 
Jo 

of untrapped particles in the system. To calculate this quantity we note that in thermal equilibrium the particle is 
equally likely to be in any cell in the system. Focusing on the sub-ensemble in which the particle is in a particular cell, 
the probability of a particle to be in the interacting region or in the noninteracting region is given by an appropriate 
integral 

/>L+2cr poo />oo />oo 

-Pint = Z^ 1 / dq dp d(f) dn p int (q,p, </>, 7r) (53) 

J L J — 00 J — 00 J — 00 

pL pOQ pOQ pOO 

P non = Z~ l I dq I dp j d<p dm p non (q,p, </>,tt) (54) 

JO J — oo J — oo J — oo 
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10 12 14 16 18 20 



FIG. 4: Diffusion constant for an ensemble of particles in equilibrium with the oscillator chain at temperature T as a function 
of the inverse temperature /3Eb, scaled by the value D® h derived in the text. Data points include representatives from all 
sixteen combinations of the two sets of parameters indicated in the figure. The solid line is the theoretical expression derived 
in the text for the behavior of the diffusion constant at low temperatures. 



over the Boltzmann factors 



Pnon (q,p,<p,n) = cxp I -p 
associated with the Hamiltonian in each region, where 



Pint (q,P, <P, 7r) = cxp (-ft y + i (tt 2 + u 2 (f) 2 ) + ot(f> ^ 



Zo 



dq 



+ 



dp 



dq 



dp 



d-K pint (q,p, </>, n) 

) 

dir p non (q,P,i>,n) 



(55) 
(56) 

(57) 
(58) 



is a normalization factor (or single cell partition function) that makes Pint + P n0 n = 1 • Note that within each region 
the integrand in (|53JI and H54JI is independent of q, so that the resulting ^-integration just gives a multiplicative factor 
proportional to the width of the corresponding region. A straightforward evaluation of the remaining integrals gives 



Pi, 



2a 



P„. 



Le 



-0E E 



Le~P E B +2a' 



(59) 



As might be expected, at high temperatures the fraction of particles in each region is simply proportional to the 
width of that region. At low temperatures, by contrast, it becomes exponentially more likely to find the particle in 
the interaction region than in the non- interacting region. In terms of these probabilities, the fraction of untrapped 
particles in the ensemble can then be written 



fu Pnt/', 



int 



Pn. 



(60) 



where /™', the fraction of those particles in the interaction region that are untrapped in equilibrium, may be obtained 
by integrating Ea. (|47|) . The result is 



j*int 

J u 



2f3E E 



1-erf K/2/3P 



£exp(-^)exp(-^ B (l-^ 



^ T(l/A,pE B ) 



de 



2f3E B ' n£ 



■ exp (—/?£) de 



|T(3/4) 



(/3P s ) 1/4 exp(-/3P B ) 



(61) 
(62) 
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which has the following limiting behavior 

f int _ . 71 
J u 



r(3/4) 

. 16 12 



(f3E B ) 1/4 cx P (-PE B ) 



15 V 7T 



,3/2 



for (3E B > 1 
for /3E B <C 1 



(63) 
(64) 



at high and low temperatures. Clearly at high temperatures the fraction of trapped particles is small, and the resulting 
diffusion constant for an ensemble of particles in thermal equilibrium will follow the behavior given in Eq.©, but 
with algebraic corrections. At low temperatures, however, Eas. l59|l . (|60|) . and l|63|l show that an exponentially small 
fraction 



fu 



rl/2 



r(3/4) 



(l3E B ) 1/4 exp (-pE B ) 



(65) 



of particles arc in untrapped states. As a result, the diffusion constant as given by Eq. I|51|) becomes exponentially 
suppressed, and will follow the form 



D ~ D 



exp(-(3E B ) 



th" 



f3E B > 1 



(66) 



where 



u t\i — —\l - • 



2a 



(67) 



In Fig. |0J) we show numerical data for the diffusion constant of an ensemble of carriers in equilibrium with an oscillator 
chain as a function of inverse temperature, plotted on a logarithmic scale to show the exponential suppression of the 
diffusion constant as T — > 0. The solid curve is the low temperature limiting form in Eq. (|66(l . 



V. CONCLUSION 



In this paper we have introduced a simple classical version of the Holstcin polaron, namely an otherwise free mobile 
particle that experiences a linear, local interaction with vibrational modes distributed throughout the medium in 
which it moves. Similar to what has been observed in quantum mechanical treatments of the problem we find that 
there are two different types of state of the combined system: self-trapped and itinerant. Self-trapped states include 
particle-oscillator states at negative energy, as well as some states at positive energies less than the polaron binding- 
energy E B . As their name suggests, self-trapped states are immobile. Itinerant polaron states, by contrast, undergo 
motion that is macroscopically diffusive. The diffusion constant for such states varies as a power of the temperature, 
with two different regimes similar to the adiabatic and non-adiabatic limits discussed in the polaron literature. At high 
temperatures the diffusion constant for untrapped particles varies with temperature as (f3E B ) ' and the transport 
mechanism is similar to what one thinks of as band transport in a solid, with long excursions before a velocity 
reversing scattering from one of the oscillators in the system. At low temperatures transport occurs via hopping of 
the carrier between different cells, with the particle spending a long time in a given cell before making a transition to 
a neighboring one. In this limit the diffusion constant for untrapped particles varies as (f3E B )~ 3 ^ 4 . 

An ensemble of particles in equilibrium with the chain exhibits a diffusion constant that is exponentially activated 
at low temperatures, with an activation energy equal to the polaron binding energy E B . The functional form, 
i.e., exponentially activated with an algebraic prefactor, is reminiscent of the one that has been derived using the 
small polaron hopping rate for a particle moving through a tight-binding band of states Q. We anticipate that any 
mechanism (not included in the present model) that would allow for a slow equilibration or exchange of itinerant and 
self-trapped particle states would not affect the functional dependence of the diffusion constant with temperature as 
we have derived in this paper. The present model we believe is simple enough that many aspects of the statistical 
mechanics of the system in equilibrium can be derived exactly, a circumstance that will facilitate an analysis of 
transport properties. In a future paper we consider noncquilibrium aspects of this model in an attempt to study, from 
a microscopic point of view, various aspects of the noncquilibrium statistical mechanics associated with it. 

Indeed, there has been continuing interest in the theoretical problem of deriving macroscopic transport laws (such 
as Ohm's law) from microscopic, Hamiltonian dynamics 0, This problem is directly related to the Hamiltonian 
description of friction and diffusion, since any such model must describe the dissipation of the energy of a particle 
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moving through a medium with which it can exchange energy and momentum, and must predict the temperature and 
parameter dependence of transport coefficients, which will depend on the model considered and on its microscopic 
dynamics. The model we have studied here is one of a class of models in which a particle, moving through a d- 
dimensional space, interacts locally with vibrational (or rotational |l4l Il5|) degrees of freedom representing the 
medium, with the total energy being conserved during the evolution of the non-linear coupled particle-environment 
system. We have investigated the emergence of microscopic chaos arising in the local dynamics of the current model 
elsewhere 18]. For another such model [ij, with an infinite number of degrees of freedom at each point in a d- 
dimensional space, it was shown rigorously that, at finite total energy, and when a constant electric field is applied, 
the particle reaches an asymptotic speed proportional to the applied field. Positive temperature states, which would 
have infinite energy, were not treated in [l|. Their rigorous analysis would be considerably more difficult. The model 
studied in the present paper is of the same type as that in , but the infinite-dimensional and continuously distributed 
harmonic systems of the previous work have been replaced by isolated, periodically-placed oscillators. This, as we 
have seen, allows the system to be studied numerically and to some extent analytically at positive temperature. 

As we have suggested in the introduction, our system can also be thought of as a ID inelastic Lorentz gas (or 
Pinball Machine). Recall that in the usual periodic Lorentz gas, circular scatterers are placed periodically on a plane, 
with the particle moving freely between them, and scattering elastically upon contact with the obstacles. For this 
system, which has been studied extensively, it has been proven that particle motion in the absence of an external 
field is diffusive 0]. The Lorentz gas, however, cannot describe the dissipation of energy of the particle into its 
environment, since all collisions are elastic. To remedy this situation, the use of a Gaussian thermostat, artificially 
keeping the particle kinetic energy constant during the motion, has been proposed ^3 an d proven to lead to Ohm's 
law 0- Rather than using a thermostat, we have chosen in the present study to impart internal degrees of freedom 
to the obstacles, and to treat explicitly the Hamiltonian dynamics of the coupled particle-obstacle system, in a spirit 
similar to Q and the more complicated two-dimensional models of [hH Il5|. This, as we have shown, leads to a rich, 
but straightforwardly characterized temperature and model parameter dependence of the diffusion constant of the 
interacting system. 
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